#args <- commandArgs(TRUE)
q1 <- 4 #eval( parse(text=args[1]) )
q0 <- 4 #eval( parse(text=args[2]) )

MakeIndex <- function(q1, q0) {
	count.stats <- 1:(q1  + q0)
	index.treat <- combn(count.stats, q1)
	index.cont <- sapply(1:ncol(index.treat), function(j) count.stats[! count.stats %in% index.treat[, j]])	
	index <- cbind(t(index.treat), t(index.cont))
	index
}

MakeRandIndex <- function(q1, q0, n.greps = 9999) {
	q <- q1 + q0
	index <- t(cbind(1:q, replicate(n.greps, sample(1:q, q))))
	index
}

OrdTest <- function(stats, k, index, q1, q0) {
	q <- q1 + q0
	M <- nrow(index)
	r.stats <- matrix(stats[index[1:M, ]], M, q)
	perm.stats <- rowMeans(r.stats[, 1:q1]) - rowMeans(r.stats[, (q1+1):q])
	orig.stat <- perm.stats[1]
	perm.stats <- sort(perm.stats)
	as.numeric(orig.stat > perm.stats[k])
}

RandTest <- function(stats, k, sims, q1, q0) {
	q <- q1 + q0
	r.stats <- cbind(stats, replicate(sims, stats[sample(1:q, q)]))
	perm.stats <- colMeans(r.stats[1:q1, ]) - colMeans(r.stats[(q1+1):q, ])
	orig.stat <- perm.stats[1]
	perm.stats <- sort(perm.stats)
	as.numeric(orig.stat > perm.stats[k])
}

sd.mat <- function(q1, q0, lc = .001) {
	q1.mat <- matrix(1, q1, q1)
	q1.mat[upper.tri(q1.mat)] <- lc
	q1.mat <- rbind(lc, q1.mat)
	q0.mat <- matrix(1, q0, q0)
	q0.mat[upper.tri(q0.mat)] <- lc
	q0.mat <- rbind(.01, q0.mat)
	cbind(q1.mat[rep(1:(q1+1), rep(q0+1, q1+1)), ], q0.mat[rep(1:(q0+1), q1+1), ])	
}

sd.mat2 <- function(q1, q0, lc = .001) {
	q1.mat <- unique(combn(c(rep(lc, q1), rep(.5, q1), rep(1, q1)), q1), MARGIN = 2)
	q1.c <- dim(q1.mat)[2]
	q0.mat <- unique(combn(c(rep(lc, q0), rep(.5, q0), rep(1, q0)), q0), MARGIN = 2)
	q0.c <- dim(q0.mat)[2]
	q1.mat <- q1.mat[, rep(1:q1.c, rep(q0.c, q1.c))]
	q0.mat <- q0.mat[, rep(1:q0.c, q1.c)]
	t(rbind(q1.mat, q0.mat))
}

grid.rev <- function(k, index, q1, q0, reps = 10000, verbose = TRUE) {
	q <- q1 + q0
	sigs <- sd.mat(q1, q0)
	rej <- rep(NA, dim(sigs)[1])
	for(j in 1:dim(sigs)[1]) {
		sig <- sigs[j, ]
		rej[j] <- mean(replicate(reps, OrdTest(rnorm(q, 0, sig), k, index, q1, q0)))
		if (verbose == TRUE & (j %% floor(dim(sigs)[1]/10)) == 0)
			cat(paste(round(100 * j/dim(sigs)[1]), "% ", sep = ""))
	}
	if (verbose == TRUE) {
		cat("", sep="\n")
		cat("Worst-case sd:", round(sigs[which.max(rej), ], 5))
		cat("", sep="\n")
		cat("Rej frequency:", max(rej))
		cat("", sep="\n")
	}
	max(rej)
}

betaa <- .1
betab <- .1
max.grid <- 40
n.greps <- 1499

q <- q1 + q0
sc <- rep(.91, q)
#sc[4] <- .96
#sc[5] <- .96
#sc[6] <- .95
#sc[7] <- .92
#sc[8] <- .91
#sc[9] <- .91
#sc[10] <- .9
#sc[11] <- .9
#sc[12] <- .9



if( choose(q, q1) < n.greps ) {
	ks <- ceiling(choose(q,q1)*sc[q1]):(choose(q,q1)-1)
	index <- MakeIndex(q1, q0)
} else {
	ks <- ceiling((n.greps + 1)*sc[q1]):n.greps
	index <- MakeRandIndex(q1, q0, n.greps)
}
if (length(ks) > max.grid) 
	ks <- ks[floor(seq(1, length(ks), length.out = max.grid))]
set.seed(48103)
res <- rep(NA, length(ks))
for(i in 1:length(ks)) {
	res[i] <- grid.rev(ks[i], index, q1, q0)
}

res <- cbind(ks, res)
output <- paste("gridstep-aws-rev-", q1, "-", q0, ".txt", sep="")
write.table(res, output)



